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On the Generalized Hermite-Based Lattice Boltzmann Construction, Lattice Sets, 
Weights, Moments, Distribution Functions and High-Order Models 

Raul Machado* 
Faculty of Engineering and the Environment, University of Southampton, Southampton, S017 IB J, United Kingdom 

The influence of the use of the generahzed Hermite polynomial on the Hermite-based lattice Boltz- 
mann (LB) construction approach, lattice sets, the thermal weights, moments and the equilibrium 
distribution function (EDF) are addressed. A new moment system is proposed. The theoretical 
possibility to obtain a high-order Hermite-based LB model capable to exactly match some first 
hydrodynamic moments thermally 1) on-Cartesian lattice, 2) with thermal weights in the EDF, 3) 
whilst the highest possible hydrodynamic moments that are exactly matched are obtained with the 
shortest on-Cartesian lattice sets with some fixed real-valued temperatures, is also analyzed. 



I. INTRODUCTION 

The lattice Boltzmann (LB) method has been used as a viable alternative for numerical simulation of (isothermal) 
fluid flovirs for more than two decades [1], [2], [3], [4], [5], [6], [7], [8]. Yet, many aspects regarding the LB method 
can be debated, such as the choice of the construction approach to build the LB model. The continuous Boltzmann 
equation can be particularly discretized in both time and phase space [9], leading to the LB equation 



pL] : Mx + c,St,t + St) = f,{x,t) + Q{f,). (1) 

u, 

• _ Eq. (1) is a discrete kinetic equation for populations fi{x, t), where i = 1, . . . , n^ and Uq is the number of discrete 

fj . lattice velocity vectors on a Cartesian grid, fi represents the probability of finding a particle with velocity c^ at 

position X and time t in lattice units [7]. Q{fi) is the collision vector. The insertion of the nonlinear Bhatnagar- 

Gross-Krook (BGK) [10] or Welander [11] colhsion model Q{f,) = -1/t(/, - f°'^) into Eq. (1) leads to the LBGK 

equation. Note that the fi inside Q{fi) is computed at time t, i.e. the LB method is explicit, r is the relaxation 

time, non-dimensionalized with St, which describes the time of a perturbed system to return to equilibrium and it 

\^ • is related to the viscosity of the fluid. LB models are usually denoted as DdQug, where d is the dimension of the 

00 ! model [6]. For the one-dimensional {d = 1) case, LB models with a lattice set of z = 1 (c.f. Fig. 1) are low- 

^^ . order, while those with z > 1 are high-order (more about this below). The space dependence is dealt by summing 

over all the nodes of the lattice. For instance, for a low-order one-dimensional LB model, the n^ = 3 and thus 

^^ ' X]r=o fi^^i' ~ ~Ci^ f2 + Cq^/o + Ci^ fi, where Cq — and M is a non-negative integer (more about this below). These 
^^ [ LB constructions with integer q values are denoted as on-Cartesian lattice models, while those with any non-integer 
Ci value are called Ojff-Cartcsian lattice (Fig. 1). The importance of the LB equations, the asymptotic convergence 
to the continuum Boltzmann-BGK equation, the comparison to the Grad 13 moment system, etc are summarized in 

[12]. 

Generally, LB modeling boils down to find an equilibrium distribution function (EDF), Z"'^, so that some hydrody- 



j^ , namic moments, e.g. Maxwell-Boltzmann (MB) (convective) M-moments 



T.r. 



"'^'"'-^<-wy&<m <^) 



are matched, c*^^^' = c---c, M-times and 6{x) is the classical exponential function, p is the density, u is the fiow 
velocity and 6 = RT, where R is the specific gas constant and T is the temperature. The right hand side of Eq. (2) 
are the MB moments from which the density, momentum density, pressure tensor, energy fiux, rate of change of the 
energy flux conservations are obtained with M = 0,1,2,3,4 respectively. Eq. (2) is well known from the literature, 
c.f. Eq. (20) in [13], and its equivalent, Eqs. (14) and (5) in [14] and [15] respectively. The link between the needed 
lattice velocities Ug to match high order moments (2) is discussed below. 

From now on, the words thermal and isothermal are usually stated in this work in conjunction with the hydrody- 
namic moments, e.g. r.h.s. of Eq. (2). By thermal means that 9 does not need to be equal to Oq in order to match 
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FIG. 1: (Color online) Schematic representation of the relationship between the Cartesian grid and the discrete lattice velocities. 

Symbols: Dot (•): Nodes in the Cartesian grid; Circle (^): Discrete lattice velocities a. a) One-dimensional Cartesian grid, 
b) One-dimensional on-Cartesian lattice, where all Ci = i, i = integer, c) One-dimensional Oj(f-Cartesian lattice, where there 
exist at least some d ^ i, i = integer and S is one of the Cartesian-lattice mismatch distance, d) One-dimensional on-Cartesian 
lattice set c = {0, ±1,±3}, i.e. it is not the shortest lattice set for D1Q5. e) One-dimensional on-Cartesian lattice set 
c = {0,±1,±2}, i.e. it is indeed the shortest lattice set for D1Q5. One-dimensional low-order LB models have z — 1, while 
one-dimensional high-order LB models have z > 1. 



(some) hydrodynamic moments, where Bq is a particular fixed value needed to m.atch certain hydrodynamic moment. 
6q = RTq. In this context, isothermal means that 9 = Oq. How this particular 6*0 should be, is presented below (e.g. 
in connection with tables III and VII). Sometimes, 6*0 is denoted as reference "temperature". Similarly, thermal and 
athcrmal (or isothermal) weights arc denoted to those weights that arc 0-dcpcndcnt and ^o-dcpcndcnt respectively. 

The essence of the main LB idea is captured by Sauro Sued in [7], [16] and strengthen in [17] with the statement: 
"Nonlinearity is local, non-locality is (a) linear; (b) exact and explicitly solvable for all time steps; (c) space discretiza- 
tion is an exact operation". Furthermore, those theoretically fulfilled conservation laws, e.g. (2), (depending on the 
chosen LB model) are mathematically matched exactly and computationally matched to machine roundoff. To the 
LBM assets can be added: inherently parallelizable, easy handling on geometries located on-Cartesian lattice, free 
of interpolations, finite difference schemes and correcting (counter) terms (i.e. with no added extra terms evaluated 
using finite-difference schemes to obtain certain desired property). 

The low-order LBGK models contain lattices suitable to reconstruct the Navier-Stokes equation close to the incom- 
pressible limit [6], [18], [19], [13]. These models fulfill the relation (2) up to M = 1 or 2, and are usually isothermals 
(i.e. 6 = 9o), when they are free of correcting counter terms. A more free 9 value can be theoretically obtained in 
these models at the expense of the existence of spurious velocity terms. These in turn can be corrected/annihilated by 
adding extra terms evaluated using finite-difference scheme (i.e. correcting counter terms). However, such approach 
does not guarantee the main LB idea. 

High-order lattice are also studied in the literature, c.f. [20], [21], [22], [23], [13], [24], [25], [26], [27], [28], [29], [30], 
[31], where some of those models are (claiming to be) capable to recover hydrodynamics beyond the Navier-Stokes 
equation. These constructions match the expression (2) up to Af-moments, M > 2, with a certain degree of accuracy. 



The last three aforecited works, [29], [30], [31], are based on the so called "entropic" lattice Boltzmann (ELB) approach 
(c.f. appendix), while the rest are what can be called Hermite-based constructions. These high-order ELB models are 
on-Cartesian lattice but isothermals LB constructions with isothermal weights and spurious velocity terms, c.f. [31] 
(more about this below). 

Some characteristics are now outlined for the aforecited Hermite-based LB models: In [20], [23], two-dimensional 
thermal models are described with discrete velocity sets but with athermal weights, c.f. tables 1 and 2 in [23]; A 
kinetic theory study is address in [21], where ofF-Cartcsian lattice sets and athermal weights are outlined in tables 1, 
2, 3 therein. A LB model with multiple relaxation time is found in [22], where the numerical verification is based on 
off-Cartesian lattice sets and athermal weights, c.f. table 1 therein; The accuracy of the (thermal) lattice Boltzmann 
is studied in [13]; A multiple relaxation time LB model is also described in [24], where three-dimensional numerical 
validations are carried out using off-Cartesian with athermal weights, c.f. table 1 therein; A finite difference scheme is 
employed in [25] in an isothermal LB model with off-Cartesian lattice with athermal weights, c.f. tables 1,2 therein. 
In [27], an on-Cartesian Hermite-based LB model is presented, but still it is based on athermal weights and a general 
construction to obtain the shortest lattice sets are not found in the literature (more about this below). Because of 
possible discontinuities at the wall, a finite difference method is chosen in [28], due to the presence of off-Cartesian 
lattice construction. In general, finite difference schemes are adopted in many high-order LB models for stability 
issues [32]. 

There exist some other alternative (high-order) LB constructions, e.g. [33], [34], [35] and subsequent works, c.f. 
[36], [37]. Unfortunately, finite difference schemes are required. Hybrid LB constructions can be added to this group. 
For instance, an LB model is proposed in [38], where mass and conservation equations are solved due to [39], whereas 
the diffusion-advection equation for the temperature is solved separately, e.g. by using finite-difference. 

It is useful to have high-order thermal LB models on-Cartesian lattice, with thermal weights (based on the final 
results that are used in the EDF), and with the shortest lattice sets when possible. Locality has been long recognized 
as an important source of efficiency in parallel computing to lower communications overhead, c.f. [40]. Therefore, for 
the sake of (parallel) computational cost, it is good to have high-order LB models with consecutive lattice sets, e.g. 
in one-dimension (Fig. 1) Ci = consecutive integers up to z, and thus with the shortest lattice sets. A computational 
cheap LB construction makes feasible to have a complete (i.e. non-reduced), or at least a less reduced lattice set, 
needed to match (some) hydrodynamic moments. The importance of weights becomes clear at walls, where the EDF 
is (almost) equal to the density-scaled weights, f^'^ = pWi{l 4- C), where C ~ function(0, u). This, due to the flow 
velocity is (almost) zero at the walls, depending of the regime (e.g. slip or non-slip flow) and C — Q when w = 0, 
regardless 9. Hence, the importance of having thermal weights is evidenced for walls with 6 ^ 6q. 

Strategies, such as (but not limited to) interpolations and/or approximations, are sometimes implemented to deal 
with this Cartesian- lattice mismatch (c.f. Fig. 1 c) ). It should be pointed out that with the use of interpolations, 
the exact matching of the conservation laws is not guarantee and/or the locality is lost for many existing schemes [7] , 
[41]. All of this to the detriment of the main LB idea. An example: In [26], the off-Cartesian lattice problem (Fig. 1) 
is tackled so that the pointwise interpolations are avoided by adopting approximations of the non-integer values to an 
appropriate (closest) lattice grid point. Their D2Q13 model with athermal weights is already (claiming to be) able to 
capture some of the microflows features, although it is recognized in [26] that a higher order LB model is definitely 
needed (to match higher order moments and to improve accuracy). Their experience uncovers that moving from 
standard D2Q9 to their approximated D2Q13 implies not much difference in the computational cost and instability, 
[42]. However, the approximation implemented in D2Q13 becomes difficult for D2Q16 and interpolations arc needed 
and thus, the computational cost increases significantly. The D2Q21 was tried, [42], with increased computational 
cost and serious instabilities despite additional interpolations. 

Because of existence of interpolations or approximations (e.g. in some previous Hermite-based LB due to off- 
Cartesian lattice models), spurious velocity terms (e.g. in ELB method due to its macroscopic description property, 
c.f. appendix) and finite difference schemes (e.g. in alternative LB models), the main LB idea is compromised in 
some of the aforementioned high-order LB models. 

In general, two main issues are addressed in this work: i) The influence of the use of the generalized Hermite 
polynomial on the Hermite-based LB construction approach, lattice sets, the thermal weights, moments and the 
equilibrium distribution function. A new moment system is proposed. This is handled in sections II and HI. ii) 
An answer is given to the following question: Is it (theoretically) possible to obtain a one-dimensional high-order 
Hermite-based LB model capable to exactly match the first hydrodynamic z-moments thermally 1) on-Cartesian 
lattice, 2) with thermal weights (based on the final results that are used in the EDF), 3) whilst the hydrodynamic 
{z -\- l)-moments are exactly matched with the shortest on-Cartesian lattice sets with some fixed real- valued 6*? This 
is handled in section HI. 

This is a theoretical work, where the necessary equations are presented in a compact yet complete form, in order 
to avoid bulky relations. Numerical studies are presented elsewhere. The approach of presenting solely theoretical 
results about LB prior numerical simulations is adopted by other authors as well, c.f. [6], [19], [43], [20], [44], [45], 
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FIG. 2: (Color online) Comparison between the classical exponential dix) and the generalized exponential function (3) with 
/i = -i2«zii. a): 61 = 1/3 ±10"^; b): Zoomed part of a). 
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II. ON THE GENERALIZED HERMITE-BASED LATTICE BOLTZMANN CONSTRUCTION 



It is always recommended to deal with a general formulation when relations are derived. From the classical MB 
moments (2), the classical exponential function &(x) is noticed. This suggest its extensions to the term e^{x), which 
is the generalized exponential function, and it is defined as 



e^(a;) = (2.t)-i/2-m wM^,/2^^{2x), 
where WM is the Whittakcr M- function [46], defined as 

WM_i/2,^(2x) = 22'^a; J(„i/2+^) (x-)r(i + /i) + 2^^^ J(i/2+^) {x)T{^ + ^i), 



(3) 



(4) 



T{n) = (n — 1)! and J^{x) is the Bessel function of the first kind, i.e. 

Mx) = Yl ^~'' " '""^' 



m=0 



LIT n^. 

m\T{m + ^ + iy^ > 



The difference between the classical exponential, 6{x), and its generalization, 6p(x), is visualized in Fig. 2 for some 
^ values. The e^(a;) < 6{x) for x > with ^ > and for x < with /x < 0. The opposite, 6p(x) > €.{x), is obtained 
for a: < with ^ > and for a; > with ^ < 0. The Gfj,{x) and thereby the 6^(x)-dependent moment system are 
reduced to their classical 6{x) and MB moment system (2) respectively when /x = 0. In this context, the generating 
function for the generalized Hermite polynomial Hn (x) is [46] 



n 



ril-'- 



The generalized Hermite polynomials Hn'{x), introduced by Gabor Szego [47], is obtained from the relations 



(5) 



(6a) 
(6b) 



where /i > —1/2 and i"(a;) is the generalized Laguerre polynomials [46]. However, the polynomials obtained from 
(6) are sometimes normalized, c.f. [48], [49], [50], [51]. The implemented normalization in this work is 

•^"^"^ = i3(M,l/2) ' (^) 

where B(x, y) = T{x)T{y)/T{x + y) is the beta function and r(a;) = {x — 1)! is the gamma function so that 

iJ^t:'(x) = Mn{l/2) Hi^\x), (8a) 

H^^+ii^) - AA„(3/2) i?(::Vi(x), (8b) 

for n > while Hq (x) = 1. The generalized Hermite polynomials used in this work are calculated from Eqs. (8). 

A. The Thermal Weights 

The generalized Hermite-based LB construction approach is proposed in this work. Based on the definition of the 
generalized Hermite polynomial, the LB construction is not valid for fi = 1/2 — n, n = 1,2,3,..., which will be 
seen when an EDF example is presented. The thermal weights are acquired so that they and the abscissas form a 
generalized Hermite quadrature. The d-dimensional weights for the LB DdQuq models are obtained from 

nq — l d 

J2w,l[H^^Ha)=A, (9) 

i=0 a 

where a = Ca,i/\/26, A = 1 for Yla ^o ('^)j i-^- generalized Hermite order ri = 0, or ^ = otherwise and a ~ {x, y, 
z } in (9). Uq is the number of discrete lattice velocity vectors. A number of n^ + 1 relations are obtained from (9), 
and the generalized Hermite order n goes from zero to nq. For simplicity, this work is focused to a one- dimensional 
(d = 1) study from now on, i.e. DlQuq. However, this is not a limitation. Two- and three-dimensional weights can 
be obtained from algebraic products of the one-dimensional weights; for instance, it is well known that the athermal 
weights Wq = 2/3 and VFj*2 = 1/6 from the one-dimensional low-order LB models can be used to construct the two 
dimensional weights Wq = W^ ■ W^ = 4/9, Wi-i = W^ ■ W^ = 1/9 and W5-8 = W^ ■ W^ = 1/36 [6]. The same 
procedure applies for the thermal weights obtained from the formulation (9) corresponding to the low- and high-order 
LB models, c.f. [15]. The result is that a = x now and terms such as Ca^i and Ua are equivalent to Cx.i and Ux or just 
simply to Ci and u. (Do not mix the z parameter seen in Fig. 1, with the axis coordinate z, which is no longer used 
in this work). The term 

^ - ^ (10) 



is now used throughout this work. The discrete lattice velocities are contained within the vector c = {—Cz, ■ ■ ■ — 
ci, 0, ci, . . . , Cz} for a d = 1 case, c.f. Fig. 1. 

The results from Eq. (9) for the DlQn^ generalized Hermite-based LB model can be formulated as 

Wo-f[(l-^], (11a) 



ri=l 



C„ 




w^^, ^ -"-^^r_^ IJ (i-4^h (lib) 

^ n— l,r 

where the Pochhammer symbol (p, -\- ^ -\- a)^, [52], is used in 

/C™ = 2™(/i+i + a)„, m = 0,1,2,..., (12) 

B™ = K,™ when a ~ and A™ ~ /C™ when a — 1, i.e. A™' = i?™+^/_B^, and implemented in the expanded Eqs. (11) 
thereafter. Note that Eq. (11a) is a polynomial in B 9 with zeros at cf, c|, . . . c^ (i.e. when B 9 = c^) and constant 
term one. For the particular case of cf = l,c\ = q,c\ = q^,.. .,c\ = q^~^ with g > 1, Eq. (11a) can be recast in 



terms of the g-Pochhammer symbols {B 6, l/q)z- A similar analysis can be done for Eq. (lib). The g-Pochhammer 
symbol is defined as 



z-l 



(a, g), = n(l -«?')' (13) 



k=0 



where (a, (?)o = 1, and reduces to the Pochhammer symbol at the limit q — >■ 1. 

By definition, the populations are non-negative. Hence, the weights are non-negative, and thereby the thermal LB 
model is valid, provided that the 9 is within a range whose extremes (and excluded) values are obtained from the 
following relations 



^(-1)^2^-^ (/i + i),_, r-'e,(c?, ci ...,cl) = 0, (14a) 

1=0 

J2i-m'-'i^^+hz-^0^-'e,{clcl,..., /, ,...,a = 0, (14b) 



excluded 



which have in turn been obtained from Eqs. (Ha) and (lib) respectively, as reformulations by means of 
ei(cf , C2, . . . , c1) and (/i -I- \)z-i, and equalized to zero. The ei(cf , c^, . . . , c^) is the ith-elementary symmetric poly- 
nomial [46], [53], and (/i + \)z-i is the Pochhammer symbol. A recurrent 9 value, obtained from Eq. (14b), is zero 
for all z > 1. The relations (11) and (14) have been algebraically computed up to D1Q13 lattice in a general form. 
For n^ > 13 values, particular cases (e.g. with ci = 1, C2 = 2 and so on) can only be tested with today's standard 
hardware and state of art of symbolic mathematics. Note that for the lattice model D1Q3, then z = 1 in Eqs. (14) 
and the theoretical range gives 9 =]0, c1/{2ij, + l)[, which can be reduced to the particular case with /i = and ci = 1, 
as it is found in the literature, c.f. [54]. Although some weights are never zero or negative for real 9 values with 
some particular lattice sets, others can become zero or negative under the same conditions. These weights are used to 
obtain the extremes values of 9 (more about the results on this part is found in section III, in connection with table 
VIII). 

B. The Equilibrium Distribution Function 

The classical Hcrmite-bascd LB construction is derived from a combination between an exponential based 
weight function and an exponential based equilibrium function [13]. The result leads to the classical EDF 
/r = W^*Er^=o^"(«)/"K^)", where a = Cc,a/V29 and b = Ua/V29 [13]. M + N < Q, where Q is the degree 
of precision of the quadrature (c.f. [27]), N > M, so that in the low-order LB model iVmax = 3, M = 2 and thus 
Q = 5, is minimum requirement of recovering the Navier-Stokes momentum equation [27]. The generating function 

for the generalized Hcrmitc polynomial Hn (x), Eq. (5), suggests the introduction of a new equilibrium distribution 
function, f^^, i.e. 

/r^.-.i:^^^%^^fey. (15) 



ri=0 \ ^ 

For the D1Q3 generalized LB model, with N ~ 2 and A^ = 3, the result is 



/r = .m|i ' ' '''■' 



1< , 1 C^.»"a 



9{2n + l) 2 9 2 6*2(2^-^1) 

1 ^n.i^n 1 Cn.iU^ \ 



2 9^i2)i+l)(2ii + 3) 2 9^(2)1 + 1) 



(16) 



where the underlined summands correspond to the extra terms due to the iV = 3. Note that the EDF (16) is not 
valid when /x = —1/2, —3/2. The thermal weights (11) for the equation model (16) are 

9{2n+l)-cl 
Wo = 2 -^ 17a 



l 9{2^i + l) 

2 c\ 



W.a = t.^-^^lt^. (17b) 
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TABLE I: (Color online) Comparison among the "entropic", classical and the /i-generalized Hermite-based one-dimensional 



lattice LB models, where fi is /f . The a is the coordinate axis. H;^-'"'' : results from the (classical or ^-generalized) Hermite- 
based construction Eq. (15) with A'^ = 2 or Af = 3 using a number of discrete lattice velocity vectors n, and /i = or Eq. (18). 
'^ili y results from the ELB construction (c.f. appendix). The matching terms to the classical Maxwell-Boltzmann (MB) 

moments are not underlined or under-wave. The single underlined terms are conditioned to Bq — Ci/3, while the under-wave 
to ci = 1. The missing MB terms are double underlined. The terms within a box are conditioned to Eq. (18). The double 
boxed summands are spurious terms. The term Qaaa matches its corresponding MB moment at low Mach number provided 
9 — 9o — Ci/3, ci = 1 and u^ ~ 0. The mass conservation (p, not shown in the table) is achieved for all the models. 



Note that the z-EDF (c.f. Eqs. (15), (16)) equals the p-scaled i-weight (Wi) when the lattice flow velocity is zero 
{ua — 0). It is easy to see that with /j, = 0, 6 = cf/3 and ci = 1, Eq. (16) and weights (17) are reduced to the classical 
Hermite-based construction of the low-order lattice Boltzmann formulations, as they are found in the literature, c.f. 
[6], [27]. See also Fig. 3, where the weights for the D1Q3 model at 6 = Oq are represented by the symbol O — 
(circle-solid). 



Model Construction 



The formulation (16), which contains a free parameter ^, is used in this section. The results from the first three 
classical MB moments, i.e. X]i=o /°'^cf^ for -^^ — Oi l?^, which corresponds to the density, momentum density and 
the pressure tensor respectively, are analyzed. 

The density J2i=o fi'^ = p is fulfilled independently of the value of 9 and p for the relation (16) and (17) with both 
N = 2 and A^ = 3. The momentum density X]i=o ft'^'^i ^ ■? ^ P'^ i^ ^^^'^ matched under the same conditions for the 



model with N 

N = 3. 



2, c.f. H 



(0),3 



and h[^}'^ 



L,^-. diiu ±±,2) ill table 1. On the other hand, the momentum density is not fulfilled when 
In the classical Hermite-based construction the issue with A'' = 3 is solved with 9 — 9q ^ 1/3, c.f. hU' in 

table I, where p = in (16) and (17). However, the difference in this work is that both the /i and 9 can be seen as 
"free parameters" . Therefore, the model can be presented with 



136'-, 



^" = '2' 



(18) 



with the condition that 9 j^ nor c^/2 so that Eqs. (18) and (16) remain valid respectively, i.e. 9 =]0,cf/2[. Note 
that with 9 = 9o = c1/3 the /x = from Eq. (18), which is the known reference "temperature" for the low-order 
(classical) lattice Boltzmann models. 

The resulting terms pUa for both with M ~ 1 and AI = 3, corresponding to the construction H^(fs , c.f. table I, 



r(0),3 



(3) 



are obtained under similar thermal 6 conditions as in HL^ , when the generalized Hermite-based LB construction is 

introduced. On the other hand, the relation (18) has no effect on the same terms for the H|2) construction. This is 
an (algebraic) improvement over the classical Hermite-based LB construction. From the results in table I for M = 2, 



H;2') and H^g?' , and the inviscid momentum flux density [55] (c.f. Eq. (5.11) therein), the pressure p = p{c\ — 26) 
is identified. The lattice "speed of sound" yields 




(19) 

The value of so called reference "temperature" 9 = 6q ~ cf /3 is required in the low-order classical Hermite-based LB 
constructions HLx' and HL^ [6]. This eliminates spurious velocity terms in their pressure tensors, which are thereby 
matched to the classical MB moment M = 2 isothermally, c.f. table I. On the other hand, the value of 9 is found in 
the /x-generalized Hermite-based LB constructions H;^-, and H;ff?' and no spurious velocity terms are seen in table 



(2) """^ ""(3) 



I. Note that with 9 = 9f) ^ cf/S, the Eq. (19) is reduced to Csound = \/cl/3. 

It is convenient to recall at this point that the physical speed of sound Csound = Vt^ and thus, a comparison with 
the classical lattice Cgound = V9 implies that 9 = 76*, i.e. 0(1 — 7) = 0, where 9 ^ so that Eq. (16) is valid, regardless 
/i. Hence, 7=1-1- 2/Din = 1, i.e. the degree of freedom of molecules Dm = 00, which is unphysical. This leads to the 
Newton's speed of sound Cgound = V9, which uses the ideal gas equation of state p = p9 found in the Euler equation 
and 9 = constant, i.e. isothermal assumption, p = p9 is found in the pressure tensor, obtained from the MB moment 
M = 2 in (2) (more about this below). On the other hand, based on (19) the result is cf — 29 = 76*, which yields 



7 + 2 SDm + 2 

Because there is no any restrictions on the 9 value for the moments M — and 1 terms, c.f. H/2) and HJgx in table 
I, then the 9 values in (20), for monoatomic molecules (Dm = 3,7 = 5/3) and diatomic molecules (Dm = 5, 7 = 7/5), 
can be 6* = 3cf/ll and 9 = 5c1/17 respectively. However, 9 = 9o = cf/S and ci = 1 are required in the term Qaaa 
, c.f. table I, so that the classical MB moment M = 3 is matched when u'^ « 0. This limitation on the Qaaa 
term is imposed on all low-order LB models found in table I, where the "entropic" (c.f. appendix), classical and the 
/i-generalized Hermite-based one-dimensional LB models are included. The strategy of having a deviation around 9q, 
as in [54], at the expense of the accuracy of the MB moment M = 3, has found no applications among practitioners 
dealing with weakly compressible flows to the best knowledge of the author. 

The results show so far that a fixed value oi 9 = 9q is required to achieve the best possible accuracy to reconstruct 
the incompressible Navier-Stokcs equation from the low-order LB models, i.e. with z = 1, in the low Mach- number 
limit when the models are free of correcting counter terms and regardless the /x value. However, this can be changed 
when z > 1, depending on the LB construction approach. For example, when z > 3 in Eqs. (11) and A'^ > 4 in Eq. 

(15) with p ^ 1/2-n, n = 1/2, 1, 2, 3, . . . the Ei=o ^ fP'^f' ^^r M = 0,1,2 and 3 gives p, j ^ pu, P = p9{l + 2p)+pu^ 
and Q = p9u{3 + 2p) + pu^ respectively. That is, a new moment system (to be denoted as M) is obtained with the 
use of the proposed /i-gcneralized Hermite-based LB construction approach, from which the classical MB moment 
system is a particular case with /x = 0. The area of application of this AA moment system will be determined by what 
is wanted to be achieved. Anyhow, it is already noticed here that based on an approach to obtain the macroscopic 
relations (e.g. method of moments) on the LBGK equation, the solution of 

dt] + d,,P + a, ( - (r - i) (^dtP + 9:.q)) = 0, (21) 

with constant 9 yields the one-dimensional Navier-Stokes equation 

dtj + d4p + pu^) - d42pj^d,u) = 0, (22) 

where p = pc^o^^^fi and z/ = (r — l/2)CgQm^jj. When the MB moments are used, Csound = V9 is obtained. On the other 
hand, Cgound — ■\/(l + 2p)9 is obtained when the aforementioned new Ai moment system is implemented. Although 
it is noted that the value oi p = 1/Dm can be extracted from 76' = (1 + 2p)9, the existence of a free p parameter can 
be useful when dealing with (on-Cartesian) lattice sizes. More about these Ai moments in section IIL 
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FIG. 3: (Color online) Weights values (i.e. populations f°'^ with u = and p = 1) and the likely shapes of their distributions 
for the one-dimensional integers lattice set c = {— c^, ■ • ■ — ci, 0, ci, . . . , Cz}. Symbols are found in table II. 



Symbols in Fig. 3 


Model 


9 


c 


O — 
(circle-solid) 


D1Q3 (with Eq. (18)) 
(z=l) 


e = cf /3 


Cl = l 


D • ■ 
(squared-dotted) 


D1Q5 with ^ = 

(z=2) 


61 = 0.5 


Cl = l 
C2=2 


A -. 

(triangle up-dashdot) 


D1Q5 with ^ = 

(z=2) 


e = o.7 


Cl = 1 
C2=2 


V - 

(triangle down-solid) 


D1Q5 with ^ = 

(z=2) 


Eq. (31) 
c.f. Eq. 40b 


Cl = 1 
C2=3 


< — 
(triangle left-solid) 


D1Q7 with ^ = 
(z=3) 


Eq. (39) 
c.f. Eq. 41 


Cl = 1 
C2,3 = 2,3 


X ■ ■ 
(cross-dotted) 


D1Q9 with ^ = 
(z=4) 


Eq. (42b) 


Cl,2 = 1,2 
C3,4 = 3, 5 


t> - - 
(triangle right-dashed) 


DIQU with ^ = 
(z=5) 


61 = 1.0 


Cl,2 = l,2 
C3,4,5 = 3, 4, 5 


V - - 

(triangle down-dashed) 


DIQU with ^ = 
(z=5) 


Eq. (43) 


Cl,2 = 1,2 
C3,4,5 =3,4,5 



TABLE II: (Color online) Symbols corresponding to Fig. 3. 
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III. ON THE HIGH-ORDER LB MODEL 

The second issue to treat in this work is whether or not there exist Herniite-based high-order lattice Bohzmann 
DlQn^ models, Uq > 5, so that they are able to fulfill the following three characteristics within a single construction: 
capable to exactly match the first hydrodynamic z-moments with free 9 values 1) on-Cartesian, 2) with thermal 
weights (based on the final results that are used in the EDF), 3) whilst the hydrodynamic (z-f- l)-moments are exactly 
matched with the shortest on-Cartesian lattice sets with some fixed 9 values. Most of the existing high-order LB 
models arc based the classical MB moment system. Because comparisons are made in this section, the final results are 
presented here first for the particular case when ^ = 0. The case when ^ 7^ is studied subsequently. Although Eq. 
(15) becomes the same relation as the one in [13] when fi is set to zero, the proposed formulation of the weights (11), 
used directly in the final EDF, are still thermal, unlike those athcrmal used/obtained in [21] (tables 1,2,3 therein), 
[22] (table 1 therein), [56] (table 1 therein), just to mention few examples. Furthermore, the lattice Ck values in the 
thermal weights (11) can be integers. 

It can be shown that the results of combining the thermal weights (11) and the relation (15) with fi — leads 
to that the MB M-moments can be thermally matched, i.e. with ^fi'^c^'i, in the DlQn^ models up to M = z. 
M = 0, 1, 2 . . . , and Uq = 3,5,7 ... . The MB (z + l)-moment is completely matched solely at a certain fixed reference 
9 = 9o value. Alternatively, the MB (z + l)-moment can be thermally matched up to the velocity term u^~^. The 
rest of the higher order 2^-moments, where Z > {z + 1), are not completely guaranteed. The implemented Hermite 
A'' = 3, 4, 5, 6, . . . order in (15) are for the one-dimensional rig = 5, 7, 9, 11, . . . respectively. The aforementioned fixed 
values of 6*0 can be obtained, e.g. for the D1Q3, D1Q5, D1Q7 and D1Q9 models, from the relation 

z 

Y^{-\f{nq - 2 . fc)!! r-'^-efc(c?, c^, . . . , c^) = 0, (23) 

fe=0 

while for the DlQll model, the term —540 is added into the 9^ part of the polynomial generated by (23), from which 
the roots are obtained. Hence, bulky expressions are avoided in the present work. efe(cf , C2, . . . , c^) is a fcth- elementary 
symmetric polynomial, k are non-negative integer numbers. 

Some examples are outlined to corroborate the aforementioned statements. The results for D1Q5 LB model are 

E /re? = P, (24a) 

i=0 

E /^'^' - Z'"' (24b) 

i=0 
n, — 1 

E/r^? = p(e+"'), (24c) 

geq^ E /^'^? = p{Qx9u + Q^u^), (24d) 

n, — 1 
^eq ^ ^ ^cq^4 ^ p^^^Ql ^ ^^Q^2 ^ ^^^4)^ (24e) 

i=0 

S^"^ E /^cf - p{Sr9''u^S^9u^ + S^vC'), (24f) 

where Eqs. (24a)-(24e) represent the density, momentum density, pressure tensor, energy flux and the rate of change 
of the energy flux respectively. The Qi, Ri and Si are the MB coefficients. Qi = 3 and Rq = 3 in Eqs. (24d) and 
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(24e) respectively. Their values for the lattiee set D1Q5 model are: 

e^jclcl^W) 3 
^3 - Z^^5 2 ^^^^ 

_ e2(c^ci,-3g) 3 
^2 - ^^^2 2' (^^^ 

e2(c^ci,-3g) 

Ai = — ^2 ' y^') 

^3 = g^5 , (28) 

i?4 = 0, (29) 

^5 = 0. (30) 

Complete Galilean invariant is achieved when Qa = 1 (c.f. the relation (24d)). Therefore, from Eq. (25) yields 



CT 



- Po — 



v/(3c^ - 'dc{f - 24cfc 
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(31) 



which is a particular case of Eq. (23) with z = 2. Hence, ci = 1 and 02 7^ 2, otherwise the reference "temperature" 
Of) (31) is complex-valued with C2 = 2. With ci = 1 and C2 = 3 in Eq. (31) yields 6*0 = 1 + \/l0/5. Although the 
MB (z + l)-moment is matched isothermally with z = 2 for the lattice set D1Q5 model, the energy flux is partially 
fulfilled thermally for low Mach number provided that u^ « can be assumed. The rest of the MB coefficients 
(26)-(28) can be obtained with 9 from Eq. (31), ci ^ I and C2 = 3, leading to i?2 = 6 = i?^^, 5i = 1 = Sf^ and 
53 = 250(7 + 2v^2~5)/(5 + V2~bf « 6.1257 ^ Sf^ = 10. Note that i?4 = ^ Rf^ = 1 and Sg = 7^ Sf^ = 1 
regardless the values of ci and C2, i.e. they are unconditioned no matching term to the MB coefficients. 

For the D1Q7 model, the results for the J2i=o fi'^c^' moments M = 0,1,2 are the same as the relations (24a)-(24c), 
while the rest are 



-)cq 



J2fPc^ - piWu + u'), (32a) 

i=0 

i?^^= J2 fPcj = pise^ + eeu^ + R^u''), (32b) 

i=0 

n<j-l 

^^q= J2 rM = p{lbe^u + S36u'^ + S^u^), (32c) 

y°q = ^ /°<icf = p{Voe^ + i/s^'w' + 14^"^ + ^6"'), (32d) 

i.e. it is complete Galilean invariant thermally. Vo = 15 in Eq. (32d). The MB coefficients for the lattice set D1Q7 
model are: 

eM.clcl ^39) + I59^e,{clcl cj) - 816^ 
^4 = ^^, , (33) 

_ e,iclclcl-3e) + 15e\i{clcl4)~459^ 
^3 - ^, , (34) 

^5 = 0, (35) 

e^jclcl cj, -36) + 15g^ei(cf , cj, cj) - 156^ 
V2 = ^, , (36) 

V4 = ^ [- 9040' -I2cl0clcl~3cl4e-3ctcl0~ 3440 (37) 

+4cl4 + cl44 - ^1^4 + 3340^4 + 3540^4 + 334d'^4 

+ lb0^4 - '^4^4 - 3c?6ic^ + 4442 + 456i4 + 150^6*2 

-900^4-^^0^ 4^ + ib4^e'^\ 

Vf, = 0. (38) 
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The rate of change of the energy flux is completely fulfilled when i?4 = 1 in Eq. (32b). Hence, a reference value of 
9 = 9q is derived from Eq. (33), which is a particular case of the relation (23) with z = 3. With ci = 1, C2 = 2, C3 = 3 
and z = 3, Eq. (23) becomes 



1 (1225 + 735V30)2/3- 245 + 70(1225 + 735V30)i 
° " 150 (1225 + 735^30)1/3 



(39) 



Complex values of ^0 appears in the D1Q7 model too, for instance, with ci = 1, C2 = 3 and C3 > 5, which can be 
avoided using combinations such as ci = 1, C2 = 2 and C3 = 3 or 4. 

The values of 6q for some lattice set cases are found in table III, which are presented as large numbers (around 
machine precision) for the sake of compassion to [31]. The accuracy of the MB coefficients conditioned to 6* = 6*0 is 
proportional to the accuracy of the 6*0 value. The MB coefficients for the lattice set DlQn, models with Uq = 5, 7, 9 
and 11 are summed up in table IV. From the D1Q5 and D1Q7 results, the u^+i velocity terms of the MB coefficients 
belonging to the MB {z + 3) moments, i.e. 6*3 and V4 respectively, are closest to their MB values when integers lattice 
Cfc are used if the ^o is obtained using a lattice velocity set c, which is as short as possible (c.f. ^o values in table 
III), . For instance, S3 w 6.1257, S3 w 5.5732 and 5*3 w 5.3532 are obtained for the D1Q5 model with ^o computed 
with the lattice sets {0,±1,±3} (c.f. (40b) in table III), {0, ±1,±4} and {0,±1,±5} respectively The last two 
lattice sets give negative weights. 14 sa 13.7497, V4 sa 10.4718, V4 w 8.8344 and V4 sa -124.9263 arc obtained for the 
D1Q7 model with 6*0 computed with the lattice sets {0,±1,±2,±3} (c.f. (41) in table III), {0,±1,±2,±(4- 5)} and 
{0, ±1,±3, ±4} respectively. The last two lattice sets give negative weights. The rest of the MB coefficients remain 
the same as they are presented in table IV when the 6*0 values in table III are implemented for these two D1Q5 and 
D1Q7 models (more about this below). 



6»o( with c = {0, ±1, ±2}) = complex valued, (40a) 

6»o( with c = {0,±1,±3}) = 1 + ^^. (40b) 

5 



do{ with c = {0, ±1, ±2 ± 3}) = 0.697 953 322 019 683 088 24. (41) 



6'o( with c = {0, ±1,±2,±3, ±4}) = complex valued, (42a) 

eo( with c = {0,±l,±2,±3, ±5}) = 0.756 080 852 594 268 582 31. (42b) 



eo( with c = {0, ±1,±2,±3,±4,±5}) = 2.123 517 542 924 955 553 8. (43) 



TABLE III: Reference values of ^ = ^0 for some one-dimensional DlQnq models and shortest on-Cartesian lattice velocity sets 
c = {0, ±ci, . . . , ±c^}, so the MB {z + l)-moment is completely fulfilled, c.f. Eq. (2.3). Eq. (40b): D1Q5; Eq. (41): D1Q7; Eq. 
(42b): D1Q9; Eq. (43): DlQll. All populations (15) with p = 1.0 and weights (11) are positive with ^ = and given 8 = 60 
values provided that: D1Q5: < |m| < 1.145; D1Q7: < |u| < 0.761; D1Q9: < |u| < 0.346; DlQll: < |«| < 1.117. 

Based on the weights, it is easy to show that the increase of z to a value larger than one leads to a heavy tail in the 
distribution. The populations f°^ with u = and p = 1 for the one-dimensional lattice DlQrig models with Uq = 3—11 
and integer c^ values are shown in Fig. 3. For easy visualization, the likely shapes of the distributions are also drawn 
by means of interpolation among the discrete weight points of the distributions. From a light-tailed distribution when 
z = 1, for the D1Q3 model (O — (circle-solid)) to heavier tails when z is progressively increased is illustrated in Fig. 
3. In addition, the increase of the "temperature" from 6 < 60 to d^ leads to smaller kurtosis as it is seen for the cases 
D1Q5 (from D • • (squarcd-dotted), A -. (triangle up-dashdot) to V — (triangle down-solid)) and DlQll (from t> 
- - (triangle right-dashed) to V - - (triangle down-dashed)). In these two cases, the peakedness of the distributions 
are significantly affected. On the other hand, the D1Q7 and D1Q9 lattice models do not show such properties when 
9 = 60, c.f. <1 — (triangle left-solid) and x • • (cross-dotted) in Fig. 3. Coincidentally, these two particular D1Q5 
and DlQll models have the largest 9o values, as seen from the outlined values of 9o « 0.33, 1.63, 0.69, 0.75, 2.12 for 
Uq — 3, 5, 7, 9, 11 respectively. That is, fatter tails are obtained with the increase of the 9 value, as observed in Fig. 
3. 
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D1Q5 


D1Q7 


D1Q9 


DlQll 


MB coeff. 


with 
Eq. (40b) 


with 
Eq. (41) 


with 
Eq. (42b) 


with 
Eq. (43) 


R2 

Si 
S3 

V2 
V4 




1 




1 
6 


1 
6 
1 
15 

10 


1 
6 
1 
15 

10 

1 

45 

15 


6 







■1 


1 
15 


15 




6.1257 


10 










1 

45 


45 




13.7497 


15 
















1 



TABLE IV: (Color online) Values of the MB coefficients for different DlQw, LB models. The MB coefficients are seen in Eqs. 
(24d), (32b)-(32d). The presented unconditioned matching terms to the MB coefficients are not underlined or single/double 
boxed. Unconditioned no matching terms are underlined. Single boxed terms are conditioned to 9 = 9o. Double boxed terms 
are conditioned to 9 = 60, but still they are no matching terms to the MB coefficients. The chosen 60 values are seen in table 

in. 

The shape of the distribution is also altered due to the flow velocity. Recall that the populations f^'^ (c.f. Eqs. 
(11) and (15) with fj, = 0) are p-scaled velocity perturbations on the weights. An example is depicted in Fig. 4 for 
the lattice set DlQll. Here, skewness is affected with the increase of the velocity. With 6 = 1.0 and maximum 
(minimum) possible lattice velocity u = 0.74 {u — —0.74), the distribution is skewed to the left (right). Negative 
populations are obtained with further increase of the lattice velocity \u\. For instance, with populations denoted as 
{/-5, /-4, /-3, /-2, /-I, /o, /i, /2, /s, /4, /s} On the lattice set and with u = 0.75 {u = -0.75) yields /_3 < (/a < 0); 
with u = 1.1 {u = -1.1) yields /_3 < and /_2 < (/a < and /2 < 0); with u = 1.35 {u = -1.35) yields 
/_3 < 0, /_2 < and /o < (/a < 0, /2 < and /o < 0); and so on. With the reference value of Eq. (43), the 
MB (z + l)-moment is now guaranteed with z = 5, leading to a new maximum (minimum) possible lattice velocity 
oi u = 1.117 (u = —1.117), and the distribution is skewed to the left (right). In this case, negative populations start 
to show up with u ~ 1.12 {u = —1.12) and it yields /_4 < (/4 < 0); with u = 1.5 {u ~ —1.5) yields /„4 < 
and /_3 < (/4 < and /a < 0); with u = 2 {u = -2) yields /_4 < 0, /_3 < and /o < (fi < 0, /a < and 
/o < 0); and so on. That is, the first negative population is obtained where the tail is longer. Hence, the first source 
of instability at a fixed 9 value shows up on the part from which the distribution is skewed to. Upon fulfilling some 
high-order hydrodynamic moments, the presence of "thermal" tails in the distributions allows capturing high velocity 
particles. 

How large the z value must be depends on the needed MB moments to be fulfilled, which in turns is determined 
by the particular case to simulate. In general, fiows with high velocities and temperature values require large z 
values. The increase of z and the (allowed) temperature values lead to longer and fatter tails, as already mentioned. 
Normally, the (asymptotic) extremes are adopted as a starting point to study models, from which their behaviors 
in between these two sides are later considered. Most of the LB works found in the literature are based on the low 
order construction, i.e. z = 1. The other extreme, z — > 00, is now considered for the present construction. The aim 
here is not to go deep into theoretical descriptions, which can derail this work from the LB method, but to have, at 
least, a general idea about some possible/expected properties of the distribution for very large z values. Then, the 
study can be possible linked to an existing theory, from which further work can be conducted elsewhere. The study 
of the tail distribution involves the use of the cumulative distribution function (CDF), -F(c„) = X]c"=-z^cfc, and 
the complementary cumulative distribution function (CCDF), F{cn) = 1 — F{cn), where c„ are integers numbers in 
the range of [— z, — (z — 1), . . . , — 1, 0, 1, . . . , (z — 1), z]. The CDF can be entirely written in terms of the off-centered 
lattice cell weights only, e.g. Eq. (lib), i.e. 



C„<-1 C„<2 

F(c„) = (-l)^i(^") J2 W,, + Hiicn) sgn{c„) ^ M^,, +Hi(c„), 



(44) 
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FIG. 4: (Color online) Populations obtained from Eqs. (11) and (15) with fi = 0, p — 1, and the likely shapes of their 
distributions for the DlQll model with integers lattice velocity set c — { — 5, —4, —3, —2, —1,0, 1, 2, 3, 4, 5}. Symbols: A — 

(triangle-solid): with 9 — 1.0 and u = 0;\> (triangle right-dashed): with 6 = 1.0 and u — 0.74; <] - - (triangle left-dashed): 

with 9 = 1.0 and u = —0.74; D — (squared-solid) : with Eq. (43) and it = 0; x - . (cross-dashdot): with Eq. (43) and u = 1.117; 
O - . (circle-dashdot): with Eq. (43) and u = -1.117. 



where Ha{cn) is the Hcavisidc step function 




(45) 



so that Ha{0) = a is also valid and sgn(c„) = 2Hi/2{cn) — 1- 

Some (zoomed) CCDF for DlQn^ models with Uq = 3,11,13,81 and 201, which correspond to z = 1,5,6,40 and 
100 respectively, are plotted in Fig. 5 for difference 9 values and /i = 0. In general, the decay to zero of the CCDF 
becomes slower when z is increased, as seen in Fig. 5, i.e. from a rather upright CCDF for D1Q3 (* - - (asterisk- 
dashed)) to a more horizontal CCDF for D1Q201 (x — (cross-solid)). Distributions with the observed characteristics 
in Fig. 5 can be long-tailed and subexponentials. 

Some basic properties of subexponcntial distributions (at infinity) are [57]: 



(1) lim E^^ = l,VO<y<oo, 

c„^oo F{Cn) 

(2) lim e'"^"F(c„) = oo,Vs>0. 



(46a) 
(46b) 



Rigorous proofs of (46) for more general functions are found in [57] and arc not repeated in this work. Here, sketched 
proofs and examples are outlined instead, as an attempt to present the subject more accessible and intuitive to LB 
practitioners. The interpretation of the relations (46) are then used as links to show some trends and properties of 
the distribution (11) throughout examples (for some finite z values). 

Recall that the present construction allows integer lattice velocities and thus y, c„ and z are considered integers, 
i.e. I < y < Cn < z. The trivial solution y = is obviously excluded. When c„ > the ffi(c„) = 1 and sgn(c„) = 1 
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FIG. 5: (Color online) Complementary cumulative distribution function (CCDF) of the weights (populations with u = 0, 
p = 1) obtained from Eqs. (11) with fi — 0, and their likely shapes for the DlQn, models with integers lattice velocity 
set c = { — 2, — (2 — 1), ... , —1, 0, 1, . . . , 2 — 1,2} and Uq = 3, 11, 13, 81 and 201, which correspond to 2 = 1, 5, 6, 40 and 100 
respectively. Symbols: H< - - (asterisk-dashed): D1Q3 with 6 — 9o — 1/3; A . . (triangle-dotted): DlQll with 6 — 1.0; D - 
. (squared-dashdot): DQll with Eq. (43); O - - (circle-dashed): DQ13 with 6 = 2.0; <> - - (diamond-dashed): D1Q81 with 
6 = 9.0; > - . (triangle right-dashdot) : D1Q201 with 6 = 12.0; x - (cross-sohd): D1Q201 with 9 = 20.0; Solid: CCDF of the 
exponential distribution. The weights are positive within machine precision. 



and then, the CCDF of distribution becomes 



F(c„) = 1 - F(c„) = { with Eq. (44) } 

/ -1 c„<z N 

= 1 - - E ^-^ + E ^-^ + 1 

-1 C„<2 

= E ^^^ - E ^^^ 

= { with W^_,, ^Wc,} 



(47a) 
(47b) 



where the result in (47), e.g. the term / which is YTc 
number of summands ns in ^^ ^^ ^-^ Wet, is so that ns 



_i_j^ VFcfc , is zero when c„ 
c„ = z. Similarly, for F{cn - v) =Yl,l 



by definition. That is, the 



-i-.W^- 



ifhen 



c„ > 0. These weights values in F(c„), (47), and in F[cn — y) when Q < y < Cn ^ Q correspond to those located 
at the extreme of the tail. These extreme weights get closer to zero when z is increased, as they arc presented (in 
the last column) in table V. Hence, for a very large c„ < z the expression (46a) becomes a ratio between zeros. 
The I'Hopital's rule can be applied to evaluate this limit. The weights (11) can be expressed in terms of elementary 
symmetric polynomials, c.f. Eqs. (14), and the derivatives of such polynomials [58] are out of the scope of this work. 
The expression (46a) is a property of slowly varying functions (at infinity). As already noted in this work, the 
observed trend in Fig. 5 is that the CCDF varies slower to zero when z is increased. For instance, from a fast varying 
CCDF for D1Q3 (* - - (asterisk-dashed)), to a slower CCDF for DlQll with 6* = 1.0 (A . . (triangle-dotted)). 
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Model 


e 


Fig. 5 


Eq. (lib) 


D1Q3 

(z=l) 


Bo = 1/3 


* - - 
(asterisk-dashed) 


W,=i = 1/6 


DIQU 

(z=5) 


1.0 


A . . 

(triangle-dotted) 


W:,-l=4 ^ 1.3 X 10"'', 

W^^5 ~ 1.6 X lO"*^ 


DIQU 

(z=5) 


Eq. (43) 


D- . 
(squared-dashdot ) 


W,-i=4 ~ 8.2 X 10"^ 
W,=5 ~ 1.6 X 10-^ 


D1Q13 

(z=6) 


2.0 


O - - 
(circle-dashed) 


W,-i=5 ~ 3.9 X 10-^ 
W,^6 ~ 5.7 X 10-"* 


D1Q81 
(z=40) 


9.0 


0-- 

(diamond-dashed) 


H^z=30-4o = a X 10\ 

a ~ one digit, 

b ~ [-25 40] 


D1Q201 
(z=100) 


12.0, 
20.0 


t> - . 

(triangle right-dashdot) 

X — 

(cross-solid) 


Vy.=go-ioo = a X 10\ 

a ~ one digit, 

b ~ [-80 100] 



TABLE V: (Color online) Values of the weights located at the extreme of the tail for DlQn, model, where n, = 3, U, 13, 81, 2001, 
c.f. Fig. 5. The symbol ~ represents "of the order of". 

which is further changed with Eq. (43) (D - . (squared-dashdot)). The variation is shown in Fig. 5 progressively, 
for D1Q13 (O - - (circle-dashed)) with 6 = 2.0, to D1Q81 with 61 = 9.0 (0 - - (diamond-dashed)), to D1Q201 with 
9 = 12.0 (I> - . (triangle right-dashdot)), which is further changed with 6 = 20.0 (x — (cross-solid)). 

Eq. (46a) can be easily rewritten as lirric^^oo ^(s Cn)/F{cn) — > 1 with s > 0. With x = e'^", Eq. (46b) becomes 
x'^ i^(ln(x)). By definition F(ln(a:)) — ?► when x — >■ c» as a result of x = e"^" and c„ — > oo. On the other hand, 
x* — >■ oo with s > for similar reasons. The exponential CDF is 



^GVipy^n ) 



1 - exp(-sc„), c„ > 0, 
0, c„ < 0. 



(48) 



The exponential CCDF with s = 1 is plotted in Fig. 5 as a solid line-curve. From the positive side of the lattice sets, 
i.e. Cj > in Fig. 5 is observed that the decay of the CCDF for the D1Q81 and D1Q201 models are slower than 
the corresponding exponential CCDF. Thereby the name of "subexponential" . Similar studies can be done for s > 0, 
s^l. ^ 

Although it can be argued that some (physical) phenomena can be described or detected by models with long-tailed 
subexponential distributions, which in turn can be linked to extreme value theory [57], further analysis is required. In 
addition, n-modal (i.e. with ri-peaks) distributions are also observed in Figs. 3 and 4. Such studies deserve separate 
works elsewhere. 

A feasible high-order LB D1Q5 model with a lattice set {co,±ci,±C2} = {0,±1,±2} has received a deserved 
attention in the literature [59], [60], [29], [31] because it would be a good model for the Navier-Stokes equation (22) 
with the shortest on-Cartesian lattice set (Fig. 1 d)). However, in previous (isothermal) models, this one-dimensional 
five velocity model {0,±1,±2}, [59], proves intrinsic unstable [60], due to complex reference "temperature" 6*0 value 
[29]. Many causes have been attributed in order to answer the reason of such instabilities. The following statement is 
found in [31]: "In some of the earlier studies, the pattern of instability of the {0, ±1, ±2} lattice, was attributed to the 
lattice Boltzmann scheme itself [61], or to the adveetion part of the LB scheme [62], [60], or to the collision of the LB 
scheme [63], or to insufficient isotropy [64]". Although the problem is identified in the afore-cited references, different 
strategies are adopted to tackle the issue, not necessarily following the on-Cartesian approach. In [29] (which follows 
the on-Cartesian approach), the blame was put on the lattice, e.g. for the {0, ±1,±2} lattice case. It is shown later 
in this work that the proposed construction in this paper is capable to have real- valued reference "temperature" 6q 
with the shortest on-Cartesian lattice sets. 

Because of the high-order LB construction in [29], [31] is reported to be limited up to DdQ(9)'' (more about this 
below), no likely trend has been described in the literature about which lattice patterns have problems when cj, 
with fc = 0, 1, 2, . . . , i.e. consecutive integers. For the moment ^ fi'^'^a ^ , a likely trend is observed in the present 
construction in which the DlQn^ models with Uq = 3, 5, 9 have complex-valued 9o when Cfc = fc, where fc = 0, 1, . . . , z, 
while the models with 71^ = 5-1-2,9-1-2 = 7, 11 and Cfc = fc have no complex- valued 6*0- The problem with the former 
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models can be avoided, for instance, by using Ck ^ k with fc = 0, 1, . . . , (2 — 1), while Cfc > (z + 1) when k ~ z. An 
example: {cq, ±ci, ±C2, ±03, ±04} = {0, ±1, ±2, ±3, ±5} for the D1Q9. These results have been algebraically tested up 
to rig = 11 in a general form. Cases with Uq > 11, under the same algebraic conditions, are computational demanding, 
which are out of the scope of this work. 

It should be pointed out that every complex-valued ^o leads to complex-valued weights and thereby to complex- 
valued populations, which is nonsense. From probability theory, populations are nonnegative real-valued. Some 
combinations of non-consecutive c^, e.g. {0, ±1,±3, ± > 5} for the D1Q7 model, can still give a complex valued 
^o- Any ^0 value, which lead to negative populations, will contribute to instability, no matter whether they are real 
valued or not. Therefore, the calculated and presented Oq values in this work (c.f. tables III and VII ) are confirmed 
to give positive weights and populations within a given flow velocity range (in lattice units). 

All the presented relations in this work have been obtained using the Hermite construction approach. Any belief 
that the aforementioned results are merely obtained from the so called "entropic" construction [31] is discarded from 
the current results. For instance, Eq. (31) is found for the same one-dimensional lattice with n^ = 5 in [65] (c.f. 
Eqs. (lO)-(ll) therein). Furthermore, the values from Eqs. (41) and (42) are the same as those obtained in [31] 
(c.f. relations (9), (C3) and (D3), (D5) therein respectively), although presented in this work with higher accuracy. 
Finally, the isothermal on-Cartcsian lattice weights found in [31] can be also obtained from the Eq. (11) (with /i = 0) 
together with their respective 60 in Eq. (23). These similarities on isothermal weights and 6^0 values should not come 
as a surprise, taking into account how the weights and Bq values can be obtained in the ELB construction to match 
MB moments (c.f. appendix A in [66]). Similar outputs can be obtained when equivalent moments are forced to be 
matched. 

Having mentioned some similarities between the actual thermal Hermite-based construction and the isothermal 
"entropic" one (reviewed) in [31], their differences are substantial and cannot be overemphasized. Unlike the Hermite- 
based construction, the ELB method relies on macroscopic equations and its physical extension beyond these descrip- 
tions is based on adding lattice velocities, at the expense of the presence of (high) powered spurious velocity terms, 
c.f. Eqs. (8), (10) and (D6) in [31] and appendix. There are no spurious velocity terms in the relations obtained from 
the current construction, c.f. (24) and (32). The one (summarized) in [31] is based on an isothermal construction, 
i.e. mass and momentum density, and the gained 6 (through constraints) is used as a manipulation tool. Hence, from 
the pressure tensor and beyond, the results are isothermal and spurious terms are obtained. Even if the Eq. (66) in 
appendix is used to construct a high-order ELB model, the procedure is still limited to mass, momentum density and 
(the trace of) the pressure tensor. Manipulations are then needed to guarantee MB moments beyond those descrip- 
tions, say by sacrificing 6, and spurious terms show up from the energy flux and beyond. The existence of spurious 
velocity terms limits any approach to w < 1 so that they can be neglected. Note that in the current construction, 
the DlQll model with fixed 9 = 9q has maximum possible velocity of m = 1.117. The approach in [31] is outlined for 
"all possible discrete velocity sets, in one dimension", and subsequently presented up to nine-velocity set. One can 
argue that higher order ELB models can be constructed, but the flow velocity has to be limited due to the existence 
of spurious velocity terms. Although the relations (11) and (14) have been algebraically obtained up to the lattice 
D1Q13 model while (15) up to the DlQll in a general form, no mathematical restrictions are imposed in this work, no 
spurious velocity terms arise for the first MB (z + l)-moments and the value of z can be theoretically larger than that. 
For instance, some weight values computed for riq > 13 are plotted in Fig. 5. There are no rational approximations 
for the u^+^ velocity term of the MB coefficient belonging to the MB (z -I- 2) moment in the current construction, 
contrary to what is found in [29] for the lattice set D1Q5 model. In the present work, these MB coefficients related 
to the u^+^ terms are zero unconditionally, c.f. i?4 = 0, 6*5 = and Ve = in table IV for D1Q5, D1Q7 and D1Q9 
models respectively. 

Note that the use of the weights (11) with (15) and fi — leads to thermally matched MB z-moments exactly, 
mathematically speaking, with integer lattice velocities Ck for the DlQrig high-order models, as already mentioned 
above. Hence, no interpolations or approximations are theoretically needed to coincide with the Cartesian grid nodes, 
c.f. Fig. 1. The present thermal construction is more general and accurate than the isothermal low Mach number 
approach (summarized) in [31]. 

Similarly, there exist a link between the computed weights from Eq. (11) with jj. ~ and those found in the 

literature by other authors, e.g. in [21]. For instance, the roots of the fifth-order Hermite polynomial iJg {x/\/2) ~ 
are ^,±\fl^~^/T^, ±yj^^sf2^. With z = 2, ci = ^5 - \/2~5 and cs = ^5 + \/2~5 in Eqs. (31) and (11) lead 
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Cocff. 


D1Q5 


D1Q7 


D1Q9 


DlQll 


Qs 
-R2 
Ri 

Si 

S3 
Ss 
V2 
Vi 




1 




1 

2(3 + 2m) 

1 


1 

2(3 + 2a<) 

1 

(3 + 2m)(5 + 2m) 

2(5 + 2m) 

1 ^ 1 


1 

2(3 + 2m) 

1 

(3 + 2m)(5 + 2m) 

2(5 + 2m) 

1 

3(3 + 2m)(5 + 2m) 

3(5 + 2m) 

1 




2(3 + 2p) 











(3 + 2m)(5 + 2m) 


(^ 


i + 2M)(5 + 2M) 




1 ■■■ 1 






2(5 + 2/.) 




















3(3 + 2m)(5 + 2m) 


3( 


3 + 2m)(5 + 2m) 




1 ■•■ 1 






3(5 + 2m) 

































TABLE VI: (Color online) Values of the M coefficients for different DlQriq LB models. The coefficients are seen in Eqs. 
(24d), (32b)-(32d). The presented unconditioned matching terms to the coefficients are not underlined or single/double boxed. 
Unconditioned no matching terms are underlined. Single boxed terms are conditioned to ^ = So, where 60 — function(ci, /^). 
Double boxed terms are conditioned to 9 = 60, but still they are no matching terms to the coefficients. Qi = (3 + 2fi) and 
i?o = (1 + 2^)(3 + 2/i) in Eqs. (24d) and (24c) respectively. Vb = (1 + 2^)(3 + 2^)(5 + 2^) in Eq. (32d). 



to 0n = 1 and 



Wo 

Wi,2 
1^3,4 




60 



(49a) 
(49b) 
(49c) 



respectively. Hence, in this context, the current construction can be reduced to the particular case of athermal 
weights in [21], where the results in Eqs. (49) are found in table 1 therein. Note that with a = v5 — \/2 • 5 and 
V5 + %/2~5 in Eq. (28) the S3 = 10 



C2 = V t) + V'2 ■ 5 in Eq. (28) the 5*3 = 10 = S'a , while the rest of the MB coefficients remain the same as they are 
presented in table IV for the D1Q5 model. That is, only the MB coefficient related to the u^+^ velocity term belonging 
to the MB (z + 3) moment is improved when the the aforementioned non-integers lattice Ck are implemented. This, 
at the price of using, for example interpolations due to the presence of Cartesian-lattice mismatches, c.f. Fig. 1. 
Thus the computational cost is increased and the main LB idea is not guaranteed. Hence, the choice of on-Cartesian 
integer lattice velocity c^ seems more appealing. 

The implementation of the classical Hermite polynomial shows that for certain lattice sets (c.f. results (40a) and 
(42a) in table HI) the hydrodynamic {z + l)-moments are not matched with the shortest on-Cartesian lattice sets 
with some fixed real-valued 9. Hence another construction is needed to accomplish it, whilst preserving the the 
on-Cartesian and non-fixed 9 value properties for the other hydrodynamic z-moments. At the end of section II, 
the pressure tensor (P) and energy flux (Q), computed from a thermal high-order ^-generalized Hermite-based LB 



construction (for DlQn 



9' 



> 7), are used to obtain the lattice pressure and kinematic viscosity. The analysis 
previously done for the MB moments is now equivalently carried out for the new M moment system. Hence, the 
advantages of the entire new proposed construction are now outlined in details, in order to answer the question in 
the second issue at the end of section I. 

Similarly to the MB coefficients, the Ai coefficients are also denoted here as Qt, Rt, St and Vi, which correspond 
to each u' velocity terms in Eqs. (24d)-(24f) and (32d) respectively. The term 9 = RT, found in the MB moment Eq. 
(24c) and in the lattice kinematic viscosity, becomes RT{1 + 2fi) in the A4 moment system (c.f. section II). That is, a 
sort of a rescaled T value, seen from an algebraic point of view. The rest of the 6'-linked coefficients, e.g. Qi, ^0, ^2, 
i?2, Si, S3, Vq, V2, V4 c.f. Eqs. (24d)-(24f) and (32d), are subjected to equivalent transformations, c.f. table VI. A 
comparison between tables IV and VI reveals that the 6'-linked MB coefficients a-9^ become b-9" ■ (c+2/.i) • ((i+2/i) • . . . , 
n-times in the corresponding A4 coefficients, where a = b ■ c- d ■ . . . and c, d, etc are positive odd integers. The other 
coefficients, e.g. Q3, R4, S5 and Vg, are the same as the MB coefficients. 

The obtained new Ai moments, whose some coefficients are found in table VI, can be directly generated in a similar 
way as the MB moments are acquired from Eq. (2). A way can be to use the following relation 



i=0 



e - 




(50) 
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from which the generated moments (rhs in (50)) are expanded and the containing terms J-^ are subsequently substi- 
tuted by 

j-mMax ^ " jT-~^ n-2m- mod{n - 2m + 1, 2) + 2^ 

-1--1- n-2m-mod(n-2TO+ 1,2) ' ^'^' 

The term mod(i, j) = k represents the modulo operation, where k is the reminder on division i/ j. Note from Eq. (51) 
that J"™"-- = 1 when ^ = and then Eq. (50) gives the MB moments (c.f. Eq. (2)). 

The equivalent procedure used to determine Eq. (23) is carried out now to obtain 9 = 9q from 



E 

A:=0 



(-1)' 2^+-^(/i + i),+i_, e^-'e,(clcl ...,cl) = 0, (52) 



for the D1Q3, D1Q5, D1Q7 and D1Q9 models, while for the DlQll model, the term -36(5 + 2^)(3 + 2fi) is added 
into the 9^ part of the polynomial generated by (52), from which the roots are obtained. The terms (/i + ^)z-\.i-k 
and efe(cf , C2, . . . , c^) are the Pochhammer symbol and the /cth-elementary symmetric polynomial respectively. When 
^ = 0, the term 2^+*~''"(^-|- ^)z+i-k and Eq. (52) become {uq — 2 • fc)!! and Eq. (23) respectively. The final results are 
presented in table VI and some examples in table VII, which become equal to those in tables IV and III respectively 
when /x = 0, and has been acquired following the same procedure around the Eqs. (25)-(31). The example values of 
^0 in table VII are presented in that way (around machine precision) so they can be compared to their corresponding 
6*0 values in table III. Note that the Oq values in table VII are lower than those in table III. The peakedness of the 
distribution is affected when ^ = is increased to some value ^ 7^ 0, c.f. Fig. 6, where the distribution of the DlQll 
model is depicted with 9 = 1.0, p = 1.0 and u = 0, for both fi = and /i = 1/10. Long tailed and subexponentials 
distributions can also be obtained for high z values when /i 7^ 0. 



15 S 

6»o( with c = {0,±1,±2} with/i = 1/3) = \- ^33. 

34 374 


(53) 


do{ with c = {0, ±1, ±2 ± 3} and ^t = 1/5) = 0.498 Oil 143 151 771 857 6. 


(54) 


eo( with c = {0, ±1,±2,±3,±4} and A* = 1/5) = 0.531822 832 492 398 970 86. 


(55) 


eo( withc = {0,±l,±2,±3,±4,±5} and At = 1/10) = 2.056 245 985 122 330 338 8. 


(56) 



TABLE VII: Examples of reference values of 60 — function(ci,/x) for some one-dimensional DlQng models (c.f. Eq. (52)) and 
shortest on-Cartesian lattice velocity sets c — {0, ±ci, . . . , ±Cz}, so the A4 {z + l)-moment is completely fulfilled. Eq. (53): 
D1Q5; Eq. (54): D1Q7; Eq. (55): D1Q9; Eqs. (56): DlQll. All populations (15) with p = 1.0 and weights (11) are positive 
with given fi and 9 ^ 60 values provided that: D1Q5: < |u| < 0.802; D1Q7: < |u| < 1.081; D1Q9: < |u| < 0.443; DlQll: 
< |k1 < 1.323. 

The results obtained from the shortest lattice (the most desirable lattice due to their "more local" property) are 
presented in table VII with /i 7^ 0. Unlike with the MB moments, the Ai moment system gives 9o = function(ci, /^), c.L 
Eqs. (23) and (52), where the extra fi parameter can give the theoretical possibility to obtain the shortest lattice with 
a real- valued 6*0. Hence, the lattice should not longer solely blamed for the existence of complex 6*9 values when 9o = 
function(ci, /^). For example, with lattice velocity integers {0, ±1, ±2} and {0, ±1, ±2, ±3, ±4}, the complex valued 6*0 
in Eqs. (40a) and (42a) become real valued in Eqs. (53) and (55) respectively with /i 7^ 0. All this while the previous 
advantages acquired from the use of the MB moments are kept, i.e. the first Ai ^-moments are thermally matched 
with the use of the /^-generalized Hermite- based LB construction on-Cartesian lattice. The M {z + l)-moment is 
isothermally fulfilled with the shortest on-Cartesian lattice set. This is clearly an advantage of the new proposed LB 
construction, compared to the previous Hermite and ELB models. 

It has already mentioned in this work that the theoretical valid range of example of 9 values, useful in the first 
hydrodynamic z-moments, can be obtained from Eqs. (14) so that the thermal weights (in the EDF) are non-negative. 
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FIG. 6: (Color online) Weights values (i.e. populations fl^"^ with u = and p = 1) and the likely shapes of their distributions 
for the DlQll with 9 — 1.0. V — (triangle down-sohd): /i = 0; [> - - (triangle right-dashed): ^ = 1/10. 



However, only the valid range of 9 for the low-order D1Q3 model is given so far (c.f. section II). LB equations are 
discrete formulations, and thus the possible values of 9 can be segmented. Only the largest ranges for each case are 
presented. The largest valid ranges of 9 for the DIQS, D1Q7, D1Q9 and DlQll lattice models are given in table VIII. 
It should be noted that the reference 6*0 values in tables III and VII arc within the theoretical valid ranges of 9 given 
in table VIII. Based on tables III and VIII, it is interesting to note that for /i = 0, the 9q value of the DlQ(5-|-n) 
lattice can become one of a extreme ^min/max value for the next DlQ(5-|-n-f2) lattice model, where n = 0,2 and 4, 
c.f. Eqs. 9o = 1/3 (in section II), (40b), (41), (42b) and (57a), (59b), (61a), (63a) respectively. Similar findings can 
be observed for some /i 7^ cases, c.f. Eqs. (54) and (62a). 

Eq. (22) is presented for 9 = constant, but the fixed 9 value is not specified. It can be argued that 6*1 = (1 + 2/i)6'2 
algebraically, for non trivial values. Although there exist many 9i, /i and ^2 values, some few concrete examples 
are now outlined. For the D1Q7 case: With 9i = 0.69721560041248060075, which is within the extremes (59) in 
table VIII, and /i = 1/5, leads to a 02 equal to the reference value (54) in table VII. For the D1Q9 case: With 
9i = 0.7445519654893585592, which is within the extremes (61) in table VIII, and ^ = 1/5, leads to a 6*2 equal to 
the reference value (55) in table VII. For the DlQll case: 9i = 1.8, which is within the extremes (63) in table VIII, 
and /.t = 1/10, leads to 9-2 = 1.5, which is within the extremes (64) in table VIII. Hence, exactly same results can be 
theoretically obtained from the Eq. (22), for both MB and A4 moment systems. The D1Q5 case can be excluded 
because it only matches the fourth (Af = 3) hydrodynamic moment (i.e. complete Galilean invariant) at 9 — 9q. 
However, with a chosen appropriate flow velocity so that m'^ « it is observed that with 9i ~ 25/34 + 5-\/33/374, 
which is within the extremes (57) in table VIII, and ^ = 1/3, leads to a 02 equal to the reference value (53) in table 
VII. This 02 = 9o value is needed to match the fourth hydrodynamic M moment, with the shortest on-Cartcsian 
lattice set. 

Although this work is predominantly theoretical a numerical test is presented, to demonstrate the feasibility of 
the proposed ^-generalized Hermite high-order LB construction for both /z = and ^ ^ 0. Two D1Q9 cases are 
chosen: a) with lattice set {0, ±1, ±2, ±3,±5}, n = 0, reference value 0o found in (42b), table HI; b) with shortest 
lattice set {0,±1,±2,±3,±4}, ^ — 1/5 and reference value 0o found in (55), table VII. With kinematic viscosity 
1/ ~ {t ~ 1/2)(1 + 2fi)9o = 1/30 and their corresponding values of fi and 9o, the r values needed in the LBGK 
formulation for each case are obtained. A one-dimensional shock tube is simulated with an initial density ratio of 1:2 
so that p = 1.0 for x < L/2, L being the leirgth of the domain, and p = 0.5 otherwise. The results are depicted at the 
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Smin( with c= {0,±1,±3} and /i = 0) = 
emax{ with c= {0,±1,±3} and A* = 0) = 


1 

3' 
= 3. 


(57a) 
{57b) 


0mm( with c = {0, ±1, ±2} and fi = 1/3) 
Smaxl with c = {0, ±1, ±2} and fi = 1/3) 


3 

IT' 
12 

ll' 


(58a) 
(58b) 


emm( with c = {0, ±1, ±2, ±3} and ^ = 0) = 
9max( with c = {0, ±1, ±2, ±3} and fi = 0) = 


5 
1+^. 


(59a) 
(59b) 


e,^in( withc = {0,±l,±2,±3} and/i = 1/5) = — - 5 "^"^ . 

0mM with c = {0, ±1, ±2, ±3} and /i = 1/5) = 1.401 283 831 980 340 563 9. 


(60a) 
(60b) 


en,in( withc= {0,±1,±2,±3,±5} and/i = 0) = 0.697 953 322 019 683 088 24. 
emax{ withc= {0,±1,±2,±3,±5} and/i = 0) = 2.881311061716 039 428 2. 


(61a) 
(61b) 


0min( withc= {0,±1,±2,±3,±4} and At = 1/5) = 0.498 011143 151771857 6. 
6»max( with c = {0, ±1, ±2, ±3, ±4} and ^l = 1/5) = 1.829 636 973 881 101 141 2. 


(62a) 
(62b) 


6lmin( withc = {0,±l,±2,±3,±4,±5} and At = 0) = 0.756 080 852 594 268 582 31. 
^maxC withc = {0,±l,±2,±3,±4,±5} and At = 0) = 2.175 382 386 573 040 694 7. 


(63a) 
(63b) 


6ln,in( withc={0,±l,±2,±3,±4,±5} and At = 1/10) = 0.963 908 781629 469 643. 
6»max( with c= {0,±1,±2,±3,±4,±5} and At = 1/10) = 2.141493 081363 463 722. 


(64a) 
(64b) 



TABLE VIII: Theoretical largest valid range of examples of 6 for the D1Q5 (Eqs. (57)-(58)), D1Q7 (Eqs. (59)-(60)), D1Q9 
(Eqs. (61)-(62)) and DlQll (Eqs. (63)-(64)) lattice models so that the weights are positive, c.f. Eqs. (14). The extremes 
should be excluded, i.e. 9 =]&min, Smax[. 

same time step in Fig. 7. In Fig. 7, a compressible front moving into the low-density region while a rarefaction front 
moving into a high-density region are observed, as expected from these kinds of simulations. The observed oscillatory 
pattern at the shock is common in the lattice Boltzmann schemes, c.f. [17], [29]. Both D1Q9 cases are on-Cartesian 
LB models, but the one with /i 7^ 0, Fig. 7 b), has the shortest lattice set {0, ±1, ±2, ±3, ±4}. The sole purpose of this 
numerical test is for a simple computational proof of concept. Further numerical studies are carried out elsewhere. 

Certainly, it would be interesting to obtain more hydrodynamic terms from the M moment system beside p = 
pRT{\ + 2fj,) and i^ ~ {t ~ 1/2)RT{1 + 2/x). However, i) it would derail this work from its two main general issues, 
mentioned in the abstract and in the last paragraphs at the end of section I, ii) such study has to be placed within 
the context of another set of references because based on previous works (c.f. [39], [67], [22], [68], [69], [24]) another 
LB formulation would be needed to compensate the limitations of the LBGK. This deserves a separate work and it is 
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FIG. 7: (Color online) Simulation of the one-dimensional shock tube problem by the ^-generalized Hermite high-order LB 
construction for two D1Q9 cases with kinematic viscosity u — 1/30: a) fi — 0, lattice set {0, ±1, ±2, ±3, ±5} and 6o found in 
(42b), table III; b) /i = 1/5, shortest on-Cartesian lattice set {0, ±1, ±2, ±3, ±4} and real-valued ^o found in (55), table VII. 
Length of the domain 8 x 10*", time step 3 x 10*", p = 3. 
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constructions 
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ii) Thermal 
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iv) Spurious 
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v) Shortest lattice 
sets 


ELB construction, c.f. [31], 
Refs. therein and appendix. 


Yes 


No 


No 


Yes 


No 


Previous Hermite 
construction, c.f. [27] 
and references therein. 


Yes 


Yes 


No 


No 


No 


Proposed 

Hermite 

construction 


with /i = 


Yes 


Yes 


Yes 


No 


No 


with /i 7^ 


Yes 


Yes 


Yes 


No 


Yes 



TABLE IX: Comparison among non-mixed high-order LB constructions (theoretically free of finite difference scheme) based 
on whether or not they can have: i) lattice velocities on-Cartesian, ii) their first hydrodynamic z-moments thermally fulfilled 
exactly, iii) thermal weights (based on the final results that are used in the EDF), it;) spurious velocity terms in their first 
hydrodynamic z-moments, v) a stable (one-dimensional) model with the shortest on-Cartesian lattice sets (e.g. Ci — consecutive 
integers, Fig. 1) capable to exactly match the hydrodynamic (z-l- l)-moments with some fixed real- valued 6. Positive properties 
are underlined. 

presented by the author elsewhere. 

Summarizing, a comparison is made in table IX among some high-order LB models. The positive properties are 
underlined. Obviously, the new proposed LB construction with ^ 7^ has the most (theoretical) advantages. It is 
noticed that the insertion of new advantageous properties are obtained whilst previous advantageous properties are 
kept. 



IV. CONCLUSION 



The /i-generalized Hermite polynomials is proposed into the lattice Boltzmann (LB) approach, where /x 7^ 1/2 — n, 

n = 1/2, 1,2,3, In the process, a new moment system (denoted as M) is proposed (c.f. Eqs. 50 and 51). The 

M moment system reduces to the Maxwell-Boltzmann (MB) moments when fj. = 0. A new equilibrium distribution 
function (EDF) based on the /j,-generalized Hermite polynomials is also introduced (c.f. Eq. (15)). The new proposed 
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higher-order LB construction is constrained into the main LB idea (c.f. [7], [16], [70], [17]). The new formulation is 
on-Cartesian lattice sets, in order to avoid the (theoretical) need of interpolations, approximations or finite difference 
schemes. A single formulation for one-dimensional thermal weights (based on the final results that are used in 
the EDF) is introduced in this work for an unlimited Uq on-Cartesian lattice grid points, where Uq = 2z + 1 and 
z = 1,2,3,... (c.f. Eqs. (11) and Fig. 1). This is in clear contrast to previous athermal weights (c.f. [21] (tables 1,2,3 
therein), [22] (table 1 therein), [56] (table 1 therein)). Two- and three-dimensional thermal weights can be obtained 
by mean of algebraic products of the one-dimensional thermal weights. The thermal term means in this work that 
the "temperature" 9 does not have to be a fixed value, unlike in the isothermal case. A fixed "temperature" 6 — 6q 
value is denoted as reference value. The EDF (c.f. Eq. (16)) is of the form f^'^ = pWi{l + C), where C is zero 
when the flow velocity (m) is zero, and the importance of the weights values Wi is noticed. The flow velocity can 
be zero at the boundaries/walls, e.g. for a laminar channel flow. The possibility of having a high-order thermal LB 
model on-Cartesian lattice with thermal weights (Wi), free of interpolations and finite different schemes, can be useful 
when dealing with boundaries, in particular for prospective models treating "heated" walls. The first hydrodynamic 
z-moments, where J^i fi^¥ ^ -^ — 0, 1, 2 . . . , z, are exactly matched thermally when the aforementioned introduced 
formulation of the thermal weights within the proposed EDF is implemented. 

Another important issue to deal with is to obtain a high-order LB construction so that it is as local as possible, and 
thus efficient (parallel) computations can be carried out. In previous higher-order LB models, some one-dimensional 
lattice sets, e.g. Uq = 5, prove intrinsic unstable when the shortest on-Cartesian lattice sets are used (c.f. [59], 
[60], [29], [31]). In the proposed high-order LB construction with ^ = 0, some complex- valued are also obtained 
in those cases (c.f. table III), as in [31]. However, since 6 ~ function(^,Ci) is obtained in the full scale proposed 
LB construction, this is changed when /U ^ (c.f. table VII and Fig. 7). The new high-order LB formulation 
proposes a general construction to obtain real-valued 9 using the shortest on-Cartesian lattice sets in one-dimension. 
Therefore, the highest hydrodynamic moments that can be exactly matched using the shortest on-Cartesian lattice 
sets in one-dimension in the proposed high-order LB construction, are the hydrodynamic {z -\- l)-moments, which are 
fulfilled isothermally (i.e. with 6 = 6q). Also, because of the thermal and accurate nature of the obtained relations, 
the presented approach is better than the one summarized in [31], where a z-limited isothermal less accurate (with 
spurious velocity terms) construction is reported (c.f. a comparison in table IX). 

A single relation (c.f. Eqs. (14)), from which valid ranges of 9 values can be extracted (c.f. table VIII for some 
ranges) so that the thermal weights are non-negative, is introduced. The reference 6q values, needed to exactly match 
the hydrodynamic (z -|- l)-moments in the proposed LB construction (for both /i = and /i 7^ 0), for some of one- 
dimensional on-Cartesian lattice sets are provided (c.f. tables III and VII). These ^o values can be obtained from a 
single relation (c.f. Eq. (52)) for the DlQng models with rig = 3 — 11, which is also introduced in this work. Valid 
ranges of fiow velocities (u) in lattice units for the DlQn^ models with n^ = 5 — 11 with some fixed 9 values, so that 
the populations are non- negative, are also presented (c.f. captions in tables III, VII and Fig. 4). It should be pointed 
out that proposed formulations for the thermal weights and the EFD (c.f. Eqs. (11) and (15)) are presented in a 
general form, where the on-Cartesian case is a particular one. A trade-off between the use of on- and Ojff- Cartesian 
lattice sets when it comes to matching hydrodynamic coefficients is described by an example. The influence of the 
temperature, z value and flow velocity on the likely shapes of the distributions are also discussed. For high-order LB 
with very high z values, the distributions can be long-tailed and subexponentials (c.f. Fig. 5). Hence, the high-order 
LB construction is put on a firm theoretical ground (c.f. [57]), from which further theoretical studies can be carried 
out. 

The asked question in the abstract, introduction and section III is answered: Yes, it is (theoretically) possible 
to obtain a high-order Hermite-based LB model able to fulfill the following three characteristics within a single 
construction: capable to exactly match the first hydrodynamic z-moments thermally 1) on-Cartesian, 2) with thermal 
weights (based on the final results that are used in the EDF), 3) whilst the hydrodynamic (z -1- l)-moments are exactly 
matched isothermally using the shortest on-Cartesian lattice sets. 
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Appendix 

For completeness, self-consistent, some of the main ideas behind construction of the ELB construction are sum- 
marized in this appendix. The derivation starts with a discrete _ff-function (Boltzmann ansatz) H = (—5) = 

X]i=o fi^^&ifi/^i)j where 5 is a concave function, which represents the entropy and 7f is a convex function so 
that H — —S [71], [43], [72]. The EDF is obtained by minimizing the H function upon the constraints (2) up to 
M ~ A/niax, where Af^ax = 1 or to A/max = 2, depending on the model, i.e. 

-Qir+ 1^ ^{M),a qJ, = 0, (65) 

where <^{M) represents the r.h.s. of Eq. (2) at M-momcnt, such that $(0) = p, $(1) = pu^ = ja and $(2) ^ pfaa = 
P(Csound + "q)- The lattice speed of sound can be Cgound = V^- xIm) a ^""^ ^^^ Af -Lagrange multipliers for the each 
of the Af = 0, 1, 2 MB-moments respectively. xIq) a ^^ X(o) since the density is a scalar quantity. Because of constants 
are immaterial [73] [74] and with X(a/),q = ^(X(m) a)' ^^^ solution for /i, which becomes the f'^'^, yields a thermal 
product form 

Q = {a;,y,z} \ \ A/=0 / / 

where d stands for the dimension of the problem. The Lagrange multipliers are found upon substituting (66) into the 
constraints (2) up to Ad = Af^ax- The result for DdQS'' can be formulated in the following form 



/r = p n ^^. 

a = {x,y, z} 




p0(Paa - Cl) 



pPaa - jaCl 



./ci 



(67) 



where the underlined term appears when Afmax = 2. Note that the insertion of 6' = ^o = cf /3, ci = 1 and ja = pua 
into the one-dimensional weights (17) with p = and also into Eq. (67) reveals the reduced formulation (3) in [75]. 
It is interesting to note that no product form is mentioned in [75]. In addition, the product form in [31] is only up to 
Afmax = 1 in Eq. (66), while no (thermal) product form is implemented in [76] (c.f. Eq. (25) therein). 

When A'/„iax = 1 in (66), the result becomes similar as in (67) but without the underlined term and fixed 9 ~ ci/S, 
Cl = 1, meaning that only the density and the momentum density are fulfilled. Therefore, there exist spurious velocity 
terms in the (trace of the) pressure tensor Paa for the ELB method with A/„iax = 1, as seen table I for A/ = 2, denoted 
by E?]^). This lack of accuracy is solved by using Afmax = 2 in (66) leading to (67), c.f. E%) i'^ table I for A/ = 2. 
However, because of the ELB method is based on macroscopic descriptions, the value of Afmax is limited, e.g. up to 
2, and thus the increase of the lattice set z > I cannot be equated with a similar increase in Afmax- This leads to 
spurious velocity terms in high-order ELB models, similar to those found in E^x but for higher order (MB) moments. 
This, even after using as a helping parameter in an effort to match (MB) moments. The exactness is lost with the 
presence of spurious terms and thus the main LB idea is not fulfilled for high-order ELB models. 
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